#ifndef FEM_DEMO_H
#define FEM_DEMO_H

/**
 * @file femdemo.h
 * @author shanyoumu
 * @brief 一个最简单有限元的实现
 * 
 *      实现的方程为possion方程 -△u = f，边界条件为 u = 0 。
 * @version 0.1
 * @date 2022-03-24
 * 
 * @copyright Copyright (c) 2022
 * 
 */

#include <iostream>
#include <vector>
#include <string>
#include <fstream>
#include <utility>
#include <eigen3/Eigen/Dense>

using std::vector;
using std::string;
using std::ifstream;
using std::ofstream ;
using std::pair;
using namespace Eigen ;

class point : public vector<double> {
public:
    point() : vector<double> (2) { };
 };

class Q1Element
{
protected:
    size_t index ;
    vector<size_t> dofs;      //  每个单元的自由度全局编号
    vector<point> vertexes;   //  每个单元都由四个顶点组成，保存的是顶点的坐标
    vector<pair<point,double> > weight; //计算积分时的权重以及点坐标，在标准单元中
public:
    Q1Element(int _acc = 2) ;
    void set_index(size_t _index) { index = _index ;} ;
    size_t get_index() { return index ;};
    const vector<size_t>& get_dofs() const { return dofs ;};
    void set_dofs( vector<size_t> _dofs) { dofs = _dofs ;};
    const vector<point> &get_vertexes() const { return vertexes ;} ;
    void set_vertexes(vector<point> _vertexes) { vertexes = _vertexes ;} ;
    void set_weight(int _acc) ; //  根据需要的精读来读取积分文件，初始化weight
    const vector<pair<point,double> > &get_weight() const { return weight ;}
/**
 ****************************************************************************************
 *  以下都是参考单元的函数
 */
    double basic_function(int _i , double _xi, double _eta);    //  第i个基函数在(_xi,_eta)点处的取值
/**
 * @brief 第i个基函数的梯度
 * 
 * @param _i 基函数的编号
 * @param _xi 坐标
 * @param _eta 坐标
 * @return pair<double,double> 第一个为偏xi，第二个为偏eta
 */
    pair<double,double> grad_basic_function(int _i, double _xi, double _eta);   
    Matrix2d get_Jacobian(point _p) ;
    Matrix2d get_inv_Jacobian(point _p) ;
/**
 * @brief 第i和j个基函数对总矩阵的贡献，位置为M(dofs[i],dofs[j])
 * 
 * @param _i 第i个基函数
 * @param _j 第j个基函数
 * @return double 
 */
    double con_Matrix_a_ij(int _i, int _j) ; 
};
class Q1_boundary_element : public Q1Element
{
private:
    vector<size_t> basic_index;
public:
    Q1_boundary_element(int _acc = 2) : Q1Element(_acc) { };
    const vector<size_t> &get_basicindex() { return basic_index;}
    void set_basicindex(vector<size_t> _index) { basic_index = _index;}
    pair<pair<size_t,size_t>,double> con_Matrix_aij(int _i, int _j);
};
/**
 * @brief 四边形网格剖分
 * 
 */

class RegularMesh
{
private:
    point pleft;
    point pright;
    size_t nx;
    size_t ny;
    vector<Q1Element> elements;
    vector<Q1_boundary_element> boundary_elements;
    vector<point> dofs;
public:
    void build_regularmesh(int _acc, point _pleft, point _pright, size_t _nx, size_t _ny);
    void fem_solution( string str, double (*f)(double,double));
    const vector<point> &get_dofs() { return dofs ;};

};


#endif 